function least_squares_intro()
% LEAST_SQUARES_INTRO 最小二乘问题介绍
% 
% 介绍最小二乘问题的基本概念、几何解释和应用背景

fprintf('=== 最小二乘问题介绍 ===\n\n');

%% 1. 最小二乘问题的定义
fprintf('1. 最小二乘问题的定义\n');
fprintf('给定超定系统 A*x = b，其中 A ∈ R^(m×n), m > n\n');
fprintf('最小二乘解是使 ||A*x - b||_2^2 最小的 x\n\n');

% 创建超定系统示例
m = 6; n = 3;
A = [1, 1, 1;
     1, 2, 4;
     1, 3, 9;
     1, 4, 16;
     1, 5, 25;
     1, 6, 36];
b = [1.1; 3.9; 9.2; 15.8; 25.1; 36.2];

fprintf('系数矩阵 A (%d×%d):\n', m, n);
disp(A);
fprintf('右端向量 b:\n');
disp(b');

fprintf('系统是超定的：%d个方程，%d个未知数\n', m, n);
fprintf('矩阵A的秩: %d\n', rank(A));

%% 2. 最小二乘解的存在性和唯一性
fprintf('\n2. 最小二乘解的存在性和唯一性\n');
fprintf('最小二乘解总是存在的\n');
fprintf('当 rank(A) = n 时，最小二乘解唯一\n\n');

% 检查A的列满秩性
if rank(A) == n
    fprintf('矩阵A列满秩，最小二乘解唯一\n');
else
    fprintf('矩阵A列不满秩，最小二乘解不唯一\n');
end

%% 3. 正规方程
fprintf('\n3. 正规方程\n');
fprintf('最小二乘解满足正规方程: A^T*A*x = A^T*b\n\n');

AtA = A' * A;
Atb = A' * b;

fprintf('A^T*A =\n');
disp(AtA);
fprintf('A^T*b =\n');
disp(Atb');

% 求解正规方程
x_normal = AtA \ Atb;
fprintf('正规方程解: x = [%.6f; %.6f; %.6f]\n', x_normal);

% 计算残差
residual = A * x_normal - b;
residual_norm = norm(residual);
fprintf('残差范数: ||A*x - b|| = %.6f\n', residual_norm);

%% 4. 几何解释
fprintf('\n4. 几何解释\n');
fprintf('最小二乘解对应于b在A的列空间上的正交投影\n\n');

% 计算投影
P = A * ((A' * A) \ A');  % 投影矩阵
b_proj = P * b;           % b在列空间上的投影

fprintf('投影矩阵的性质:\n');
fprintf('P^2 - P 的范数: %e (应该为0)\n', norm(P^2 - P, 'fro'));
fprintf('P^T - P 的范数: %e (应该为0)\n', norm(P' - P, 'fro'));

% 验证正交性
orthogonal_error = (b - b_proj)' * A;
fprintf('残差与列空间的正交性: ||r^T*A|| = %e\n', norm(orthogonal_error));

%% 5. 条件数的影响
fprintf('\n5. 条件数的影响\n');
fprintf('A^T*A的条件数是A的条件数的平方\n\n');

cond_A = cond(A);
cond_AtA = cond(AtA);

fprintf('cond(A) = %.2e\n', cond_A);
fprintf('cond(A^T*A) = %.2e\n', cond_AtA);
fprintf('cond(A)^2 = %.2e\n', cond_A^2);

% 演示条件数恶化的影响
fprintf('\n条件数恶化的影响:\n');
delta_b = 1e-6 * randn(size(b));  % 小扰动
b_perturbed = b + delta_b;

x_perturbed = (A' * A) \ (A' * b_perturbed);
rel_error_b = norm(delta_b) / norm(b);
rel_error_x = norm(x_perturbed - x_normal) / norm(x_normal);

fprintf('右端相对扰动: %e\n', rel_error_b);
fprintf('解的相对误差: %e\n', rel_error_x);
fprintf('误差放大因子: %.2e\n', rel_error_x / rel_error_b);
fprintf('理论上界: %.2e\n', cond_AtA);

%% 6. 多项式拟合示例
fprintf('\n6. 多项式拟合示例\n');
fprintf('用2次多项式拟合数据点\n\n');

% 生成数据点
x_data = linspace(0, 5, 10)';
y_true = 2 + 3*x_data - 0.5*x_data.^2;  % 真实函数
noise = 0.5 * randn(size(x_data));      % 噪声
y_data = y_true + noise;                 % 观测数据

fprintf('数据点数: %d\n', length(x_data));
fprintf('多项式次数: 2\n');

% 构造Vandermonde矩阵
V = [ones(size(x_data)), x_data, x_data.^2];
fprintf('Vandermonde矩阵条件数: %.2e\n', cond(V));

% 最小二乘拟合
coeffs = V \ y_data;
fprintf('拟合系数: [%.6f, %.6f, %.6f]\n', coeffs);
fprintf('真实系数: [2.000000, 3.000000, -0.500000]\n');

% 计算拟合质量
y_fit = V * coeffs;
rms_error = sqrt(mean((y_data - y_fit).^2));
fprintf('均方根误差: %.6f\n', rms_error);

% 绘制结果
figure('Name', '多项式拟合');
x_plot = linspace(0, 5, 100)';
V_plot = [ones(size(x_plot)), x_plot, x_plot.^2];
y_plot = V_plot * coeffs;

plot(x_data, y_data, 'ro', 'MarkerSize', 8, 'LineWidth', 2);
hold on;
plot(x_plot, y_plot, 'b-', 'LineWidth', 2);
plot(x_plot, 2 + 3*x_plot - 0.5*x_plot.^2, 'g--', 'LineWidth', 2);
xlabel('x');
ylabel('y');
title('最小二乘多项式拟合');
legend('观测数据', '拟合曲线', '真实曲线', 'Location', 'best');
grid on;

%% 7. 线性回归示例
fprintf('\n7. 线性回归示例\n');
fprintf('多元线性回归: y = β₀ + β₁x₁ + β₂x₂ + ε\n\n');

% 生成回归数据
n_samples = 50;
x1 = randn(n_samples, 1);
x2 = randn(n_samples, 1);
beta_true = [2; 1.5; -0.8];  % 真实参数
epsilon = 0.3 * randn(n_samples, 1);  % 误差项

X = [ones(n_samples, 1), x1, x2];  % 设计矩阵
y = X * beta_true + epsilon;

fprintf('样本数: %d\n', n_samples);
fprintf('特征数: %d (包括截距)\n', size(X, 2));

% 最小二乘估计
beta_hat = X \ y;
fprintf('估计参数: [%.6f; %.6f; %.6f]\n', beta_hat);
fprintf('真实参数: [%.6f; %.6f; %.6f]\n', beta_true);

% 模型评估
y_pred = X * beta_hat;
residuals = y - y_pred;
R_squared = 1 - sum(residuals.^2) / sum((y - mean(y)).^2);

fprintf('R²: %.6f\n', R_squared);
fprintf('残差标准差: %.6f\n', std(residuals));

%% 8. 病态最小二乘问题
fprintf('\n8. 病态最小二乘问题\n');
fprintf('当A接近列相关时，最小二乘问题变得病态\n\n');

% 构造病态矩阵
A_ill = [1, 1, 1.0001;
         2, 2, 2.0001;
         3, 3, 3.0001;
         4, 4, 4.0001;
         5, 5, 5.0001];
b_ill = [1; 2; 3; 4; 5];

fprintf('病态矩阵 A:\n');
disp(A_ill);
fprintf('条件数: %.2e\n', cond(A_ill));

% 求解
x_ill = A_ill \ b_ill;
fprintf('最小二乘解: [%.6f; %.6f; %.6f]\n', x_ill);

% 添加小扰动
delta_b_ill = 1e-6 * [1; -1; 1; -1; 1];
x_ill_pert = (A_ill' * A_ill) \ (A_ill' * (b_ill + delta_b_ill));

rel_change = norm(x_ill_pert - x_ill) / norm(x_ill);
fprintf('小扰动导致的解的相对变化: %.2e\n', rel_change);

%% 9. 加权最小二乘
fprintf('\n9. 加权最小二乘\n');
fprintf('当观测数据有不同的可靠性时，使用加权最小二乘\n\n');

% 权重矩阵（对角矩阵）
weights = [1; 2; 1; 3; 1; 2];  % 不同观测的权重
W = diag(weights);

fprintf('权重向量: [%g, %g, %g, %g, %g, %g]\n', weights);

% 加权最小二乘解
% min ||W^(1/2)(Ax - b)||²
A_weighted = sqrt(W) * A;
b_weighted = sqrt(W) * b;

x_weighted = A_weighted \ b_weighted;
fprintf('加权最小二乘解: [%.6f; %.6f; %.6f]\n', x_weighted);
fprintf('普通最小二乘解: [%.6f; %.6f; %.6f]\n', x_normal);

% 比较残差
residual_weighted = A * x_weighted - b;
weighted_residual_norm = sqrt(residual_weighted' * W * residual_weighted);

fprintf('加权残差范数: %.6f\n', weighted_residual_norm);
fprintf('普通残差范数: %.6f\n', residual_norm);

end